import numeric from "numeric"

/**
 * 粒子群算法的javascript实现
 * @author 小熊
 * @date 2018-04-27
 */

function IQPSO(dimension, ub, lb, fitness_function, particle_num, maxStep, fit_func_caller = null) {
    // 粒子个数
    this.particle_num = particle_num;
    // 单个粒子的维度
    this.dimension = dimension;
    // 单个粒子的每个维度的上限
    this.ub = ub;
    // 单个粒子的每个维度的下限
    this.lb = lb;
    // 单个粒子的适应度函数
    this.fitness_func = fitness_function;
    // 适应度函数的调用者
    this.fit_func_caller = fit_func_caller;
    // 粒子群算法中的粒子列表
    this.particle_list = [];
    // 全局粒子的最优位置
    this.gbest = null;
    // 全局粒子最优的适应度值
    this.gbestValue = Infinity;
    // 当前步数
    this.currentStep = 0;

    this.maxStep = maxStep;

    this.mbest = [];

    // 当前所有粒子的适应度值数组
    this.f = [];

    /**
     * 粒子群整体初始化
     */
    this.init = function () {
        console.log(`种群数量为:${this.particle_num}`)
        // 按照最大迭代次数开始更新
        for (let i = 0; i < this.particle_num; i++) {
            let particle = new Particle(this.dimension, this.ub, this.lb, this.fitness_func, this.fit_func_caller);
            particle.init();
            this.particle_list.push(particle);
        }
        this.update_gbest();

        this.mbest=this.calcMbest();

    };

    this.calcMbest =function() {
        //let sumbest = [0,0,0,0,0,0,0,0];
        let sumbest = numeric.rep([this.dimension],0);
        for (let particle of this.particle_list) {

            sumbest = numeric.add(sumbest, particle.pbest);

        }
        return numeric.mul(sumbest,1/this.particle_num);
    }


    this.calcInertiaWeight = function() {
        // 惯性权重的两个值
        const a0 = 1.0;
        const a1 = 0.5;
        // 根据当前迭代数和最大迭代数得到当前惯性权重
        const currentA = a0 - ((a0 - a1) * this.currentStep / this.maxStep);
        //const currentA = 0.5 * (1.5 + 0.5 * Math.cos(this.currentStep / this.maxStep) );
        return currentA
    }


    /**
     * 粒子群按照指定的步数更新
     * @param steps 整型,指定的步数
     * @returns {*[]} 列表,长度为2,第0位是全局最优位置,第1位是全局最优适应度
     */
    this.step_run = function (steps) {
        // 首先判断是否第一次运行,是的话就初始化粒子群
        //let mbest = numeric.rep([this.dimension],0);
        if (this.particle_list.length === 0) {
            this.init();
        }
        //mbest = this.calcMbest();

        // 按照指定的步数更新PSO算法
        for (let i = 0; i < steps; i++) {
            this.currentStep++;
            // 计算线性系数
            const a = this.calcInertiaWeight();
            // 更新每个粒子
            for (let particle of this.particle_list) {
                particle.update(a, this.gbest, this.mbest);

                this.f.push(particle.fitnessFunc(particle.position));
            }

            ////是否混沌扰动
            //let sigma = this.varianceArr(this.f);
            //if (sigma<0.0035){
            //    this.chaos();
            //}

            // 更新全局最优
            this.update_gbest();

            this.mbest = this.calcMbest();
        }


//打印
        console.log(`迭代${this.currentStep}次, 当前全局最优值为${this.gbestValue}`);
        return {
            currentStep: this.currentStep,
            gbest: this.gbest,
            gbestValue: this.gbestValue
        };

        // return [this.gbest,this.gbestValue];

    };


    this.varianceArr = function (arr) {
        let s,
            ave,
            sum = 0,
            sums = 0,
            len = arr.length;
        for (let i = 0; i < len; i++) {
            sum += Number(arr[i]);
        }
        ave = sum / len;
        for(let i = 0; i < len; i++){
            sums+=(Number(arr[i])- ave)*(Number(arr[i])- ave);
        }
        s=(sums/len).toFixed(4);
        return s;
    }

    this.chaos = function () {
        let randomId = Math.floor(Math.random() * this.particleList.length);
        this.particle_list[randomId].change();
        //for (let particle of this.particle_list) {
        //    particle.change();
        //}
    }


    // 更新全局最优
    this.update_gbest = function () {
        // 遍历所有的粒子历史最优解
        //let sumbest = numeric.floor(numeric.random([this.dimension]));
        //let sumbest = [0,0,0,0,0,0,0,0];
        for (let particle of this.particle_list) {
            //let pbest = particle.pbest;
            //sumbest = numeric.add(sumbest, pbest);
            if (particle.pbestValue < this.gbestValue) {
                this.gbestValue = particle.pbestValue;
                this.gbest = particle.pbest;
            }
        }


    }
}


function Particle(dimension, ub, lb, fitness_function, fit_func_caller) {
    // 单个粒子的维度
    this.dimension = dimension;
    // 单个粒子的每个维度的上限
    this.ub = ub;
    // 单个粒子的每个维度的下限
    this.lb = lb;
    // 单个粒子的适应度函数
    this.fitnessFunc = fitness_function;
    // 适应度函数的调用者
    this.fitFuncCaller = fit_func_caller;
    // 当前粒子的位置
    this.position = [];
    // 当前粒子的速度
    //this.velocity = [];
    // 当前随机点位置
    this.p = [];
    // 当前随机高斯点位置
    this.g = [];
    // 单个粒子的最优位置
    this.pbest = null;
    // 单个粒子的最优适应度值，初始值为正无穷大
    this.pbestValue = Infinity;


    // 粒子初始化公式
    this.init = function () {
        // 随机初始化位置
        this.position = this.get_rand_vector();
        // 随机初始化速度,全为0
        //this.velocity = numeric.mul(numeric.random([this.dimension]), 0);
        // 当前粒子的pbest就是这个位置
        this.pbest = this.position;
        // 当前粒子的best适应度值就是该位置的评价值
        this.pbestValue = this.evaluate();
    };

    // 获得在lb和ub范围内的随机向量
    this.get_rand_vector = function () {
        // 先产生一组0-1随机值
        //let rand_vec = numeric.random([this.dimension]);

        //fixme 混沌初始化？
        let rand_vec_0 = numeric.random([this.dimension]);
        let rand_vec = numeric.mul(4,rand_vec_0,numeric.sub(1,rand_vec_0));

        // 下面完成的过程是： lb + (ub - lb) * 初始值
        let range = numeric.sub(this.ub, this.lb);
        rand_vec = numeric.add(this.lb, numeric.mul(rand_vec, range));
        return rand_vec;
    };


    //chaos
    //this.change = function () {
    //    // 先产生一组0-1随机值
    //    //let rand_vec = numeric.random([this.dimension]);
    //    //let rand_vec_0 = numeric.random([this.dimension]);
    //    let pbest = numeric.mul(this.pbest, 1/(numeric.sub(ub, lb)));
    //    //let pbestz = numeric.mul(numeric.sub(pbest,rand_vec),1/numeric.sub(rand_vec_0,rand_vec));
    //    let pbestz = [];
    //    let pbestzValue = [];
    //    pbestz[0]= pbest;
    //    pbestzValue[0]= this.fitnessFunc(pbestz[0]);
    //    for(let i = 1; i < 30; i++){
    //        pbestz.push(numeric.mul(4,pbestz[i-1],numeric.sub(1,pbest)));
    //        pbestzValue.push(this.fitnessFunc(pbestz[i]));
    //    }
    //    const y = pbestzValue.sort((a, b) => a - b)[0];
    //    if (y<this.fitnessFunc(pbest)){
    //        pbest=y;
    //    }
    //    return numeric.mul(pbest, (numeric.sub(ub, lb)));
    //};


    // 粒子更新公式
    this.update = function (a, gbest_position, mbest) {
        // 粒子更新为更新位置
        // 更新速度：v(t+1) = 0.9*v(t) + 2*rand*(pbest-x) + 2*rand*(gbest-x)
        //fixme
        //let sumbest=[0,0,0,0,0,0,0,0];
        //for (let particle of this.particle_list) {
        //    sumbest=numeric.add(sumbest,particle.pbest);
        //}
        //const mbest=numeric.mul(sumbest,1/this.particle_num);

        let fi1 = numeric.random([this.dimension]);
        let fi2 = numeric.random([this.dimension]);
        //let fi2 = numeric.sub(1,fi1);
        //let pbest_part = numeric.mul(fi, this.pbest);
        //let gbest_part = numeric.mul(1-fi, this.gbest);
        let p_fenzi = numeric.add(numeric.mul(2, fi1, this.pbest),numeric.mul(2.1,fi2, gbest_position));
        let p_fenmu = numeric.add(numeric.mul(2, fi1),numeric.mul(2.1, fi2));
        this.p = numeric.div(p_fenzi, p_fenmu);
        //for(let i = 0; i < this.dimension; i++){
        //    this.g[i] = getNumberInNormalDistribution(this.p[i], mbest[i]-this.pbest[i]);
        //}

        this.g = getNumberInNormalDistribution(this.p, numeric.sub(mbest,this.pbest));

        //高斯变异
        function getNumberInNormalDistribution(mean,std_dev){
            return numeric.add(mean,(numeric.mul(randomNormalDistribution(),std_dev)));
        }
        function randomNormalDistribution(){
            let u=0.0, v=0.0, w=0.0, c=0.0;
            do{
                //获得两个（-1,1）的独立随机变量
                u=Math.random()*2-1.0;
                v=Math.random()*2-1.0;
                w=u*u+v*v;
            }while(w==0.0||w>=1.0)
            //这里就是 Box-Muller转换
            c=Math.sqrt((-2*Math.log(w))/w);
            //返回2个标准正态分布的随机数，封装进一个数组[u*c,v*c]返回
            //当然，因为这个函数运行较快，也可以扔掉一个
            //进一步的柯西分布
            //return Math.tan(Math.PI*(u*c-0.5));
            return u*c;
        }


        let u=numeric.random([this.dimension]);
        let v=numeric.mul(-1,numeric.log(u));
        let b=numeric.mul(a,numeric.abs(numeric.sub(mbest,this.position)));

        let r=numeric.pow(-1,numeric.ceil(numeric.add(0.5,numeric.random([this.dimension]))));
        //this.position = numeric.add(this.p, numeric.mul(r,b,v));
        //fixme
        this.position = numeric.add(this.g, numeric.mul(r,b,v));

        // 如果更新后的位置大于ub，则放到ub上
        let gt_parts = numeric.gt(this.position, ub);
        for (let i = 0; i < gt_parts.length; i++) {
            if (gt_parts[i]) {
                this.position[i] = this.ub[i] * (0.7 + 0.3*Math.random());
                //this.position[i] = this.ub[i] ;
            }
        }
        // 如果更新后的位置小于lb，则放到lb上
        let lt_parts = numeric.lt(this.position, lb);
        for (let i = 0; i < lt_parts.length; i++) {
            if (lt_parts[i]) {
                this.position[i] =  this.lb[i] + 0.3*Math.random() * (this.ub[i] - this.lb[i]);
                //this.position[i] =  this.lb[i] + Math.random() * (this.ub[i] - this.lb[i]);
                //this.position[i] =  this.lb[i];
            }
        }

        // 如果更新后的适应度值优于pbest值，则更新粒子历史最优
        let fitness_value = this.evaluate();
        if (fitness_value < this.pbestValue) {

            this.pbest = this.position;
            this.pbestValue = fitness_value;
        }
    };

    // 粒子评价函数
    this.evaluate = function () {
        // 如果评价函数的参数只需要一个,则直接传入位置
        if (this.fitFuncCaller === null) {
            return this.fitnessFunc(this.position);
        } else {
            return this.fitnessFunc.call(fit_func_caller, this.position);
        }
        // 直接调用传入的适应度函数

    }
}

export default IQPSO;